function loadaggdata(SampleStart,SampleEnd)

    dataDir = "$(pwd())/data/"

    GDPpc_data       = readdlm(dataDir * "GDPpc.csv", ',', Float64, '\n'; skipstart = 1);
    TFPgr_data       = readdlm(dataDir * "TFP.csv", ',', Float64, '\n'; skipstart = 0);
    unrate_data      = readdlm(dataDir * "UNRATE_CPS_FRED.csv", ',', Float64, '\n'; skipstart = 1);

    # initial transformations
    period_GDP = GDPpc_data[:,1]
    GDPpc      = log.(GDPpc_data[:,2])

    period_TFP = TFPgr_data[:,1]
    TFPgr      = TFPgr_data[:,2]

    period_UNR = unrate_data[:,1]
    UNR        = unrate_data[:,2]

    # compute GDP growth rates
    GDPpcgr        = zeros(size(GDPpc))
    GDPpcgr[2:end] = 400*(GDPpc[2:end]-GDPpc[1:end-1])

    # accumulate TFP growth rates
    TFP            = cumsum(TFPgr/400,dims=1)

    # trim the sample
    period_TFP_ind = (SampleStart .<= period_TFP .<= SampleEnd)
    period_GDP_ind = (SampleStart .<= period_GDP .<= SampleEnd)
    period_UNR_ind = (SampleStart .<= period_UNR .<= SampleEnd)

    GDPpc   = GDPpc[period_GDP_ind]
    GDPpcgr = GDPpcgr[period_GDP_ind]
    TFP     = TFP[period_TFP_ind]
    TFPgr   = TFPgr[period_TFP_ind]
    unrate  = UNR[period_UNR_ind];

    period_agg = period_UNR[period_UNR_ind]

    TFPgrdev     = TFPgr .- mean(TFPgr,dims=1)
    GDPpcgrdev   = GDPpcgr .- mean(GDPpcgr,dims=1)
    UNRdev       = unrate .- mean(unrate,dims=1);
    agg_data     = [TFPgrdev GDPpcgrdev UNRdev]
    n_agg        = size(agg_data)[2]

    mean_unrate = mean(unrate,dims=1)[1]

    return agg_data, period_agg, mean_unrate

end

##

function loaddensdata(SampleStart,SampleEnd,K,nfVARSpec)

    sNameFile = "K" * string(K) * "_fVAR" * nfVARSpec
    PhatDensCoef_factor = readdlm(loaddir * sNameFile * "_PhatDensCoef_factor.csv", ',', Float64, '\n'; skipstart = 1);
    PhatDensCoef_lambda = readdlm(loaddir * sNameFile * "_PhatDensCoef_lambda.csv", ',', Float64, '\n'; skipstart = 1);
    PhatDensCoef_mean   = readdlm(loaddir * sNameFile * "_PhatDensCoef_mean.csv", ',', Float64, '\n'; skipstart = 1);
    PhatDensCoef_mean_allt = readdlm(loaddir * sNameFile * "_PhatDensCoef_mean_allt.csv", ',', Float64, '\n'; skipstart = 1);
    period_Dens         = readdlm(loaddir * sNameFile * "_DensityPeriod.csv", ',', Float64, '\n'; skipstart = 1);
    MDD_GoF             = readdlm(loaddir * sNameFile * "_MDD_GoF.csv", ',', Float64, '\n'; skipstart = 1);
    Vinv_all            = readdlm(loaddir * sNameFile * "_Vinv_all.csv", ',', Float64, '\n'; skipstart = 1);
    N_all               = readdlm(loaddir * sNameFile * "_N_all.csv", ',', Float64, '\n'; skipstart = 1);

    period_Dens_ind     = dropdims((SampleStart .<= period_Dens .<= SampleEnd),dims=2)
    period_Dens         = period_Dens[period_Dens_ind]
    PhatDensCoef_factor = PhatDensCoef_factor[period_Dens_ind,:]
    PhatDensCoef_mean_allt = PhatDensCoef_mean_allt[period_Dens_ind,:]
    MDD_GoF             = MDD_GoF[period_Dens_ind]
    Ktilde              = size(PhatDensCoef_lambda)[1]

    # Load ME covariance matrices
    # Needs to be adjusted for Lambda
    VinvLam_all = zeros(Ktilde, Ktilde, sum(period_Dens_ind))
    tt_sel = 1
    for tt = 1:size(period_Dens_ind)[1]
        if period_Dens_ind[tt] == true # restrict to periods between SampleStart and SampleEnd
            Vinv_t = Symmetric(Vinv_all[K*(tt-1)+1:K*tt,:])
            VinvLam_all[:,:,tt_sel] = PhatDensCoef_lambda*Vinv_t*PhatDensCoef_lambda'*N_all[tt]
            tt_sel = tt_sel + 1
        end
    end

    return PhatDensCoef_factor, MDD_GoF, VinvLam_all, period_Dens_ind, PhatDensCoef_lambda, PhatDensCoef_mean, PhatDensCoef_mean_allt

end
